﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SmartMathLibrary
{
    /// <summary>
    /// This class provides special functions which are well known in the mathematic analysis.
    /// </summary>
    [Serializable]
    public static class SpecialFunctions
    {
        /// <summary>
        /// The standard integration value for the special functions.
        /// </summary>
        private const int StandardIntegrationValue = 1000000;

        /// <summary>
        /// Approximation of the error function.
        /// </summary>
        /// <param name="x">The x value for evaluation.</param>
        /// <returns>The approximated value of the error function</returns>
        public static double ErrorFunctionApproximation(double x)
        {
            const double a = -((8*(Math.PI - 3))/(3*Math.PI*(Math.PI - 4)));

            double numerator = (4/Math.PI) + (a*Math.Pow(x, 2));
            double denominator = 1 + (a*Math.Pow(x, 2));

            return Math.Sqrt((1 - Math.Exp(-Math.Pow(x, 2)*(numerator/denominator))));
        }

        /// <summary>
        /// The Runge function.
        /// </summary>
        /// <param name="x">The specified x position for the function.</param>
        /// <returns>The value f(x) at the specified position x.</returns>
        public static double RungeFunction(double x)
        {
            return 1/(Math.Pow(5*x, 2));
        }

        /// <summary>
        /// The Heaviside function.
        /// </summary>
        /// <param name="x">The specified x position for the function.</param>
        /// <returns>The value f(x) at the specified position x.</returns>
        public static int HeavisideFunction(double x)
        {
            if (x <= 0)
            {
                return 0;
            }

            return 1;
        }

        /// <summary>
        /// The Gudermann function.
        /// </summary>
        /// <param name="x">The specified x position for the function.</param>
        /// <returns>The value f(x) at the specified position x.</returns>
        public static double GudermannFunction(double x)
        {
            return 2*Math.Atan(Math.Tanh(x/2));
        }

        /// <summary>
        /// The Lorentzcurve function.
        /// </summary>
        /// <param name="x">The x value for the function.</param>
        /// <returns>The value f(x) at the specified position x.</returns>
        public static double LorentzCurveFunction(double x)
        {
            return (1/(Math.PI*(1 + Math.Pow(x, 2))));
        }

        /// <summary>
        /// The Lorentzcurve function.
        /// </summary>
        /// <param name="x">The x value for the function.</param>
        /// <param name="t">The centrum t.</param>
        /// <param name="s">The breadthsparameter s.</param>
        /// <returns>The value f(x) at the specified position x.</returns>
        public static double LorentzCurveFunction(double x, double t, double s)
        {
            return (1/Math.PI)*(s/(Math.Pow(s, 2) + (Math.Pow(x - t, 2))));
        }

        /// <summary>
        /// The Neilscheparabola function.
        /// </summary>
        /// <param name="x">The x value for the function.</param>
        /// <returns>The value f(x) at the specified position x.</returns>
        public static double NeilscheParabola(double x)
        {
            return (x*Math.Sqrt(x));
        }

        /// <summary>
        /// The Normal distribution function (Gauss).
        /// </summary>
        /// <param name="x">The x value for the function.</param>
        /// <param name="o">The standard degression o.</param>
        /// <param name="m">The expected value m.</param>
        /// <returns>The value f(x) at the specified position x.</returns>
        public static double NormalDistributionFunction(double x, double o, double m)
        {
            return Math.Pow((1/(o*Math.Sqrt(2*Math.PI))), (-0.5*Math.Pow(((x - m)/o), 2)));
        }

        /// <summary>
        /// The density function of the standard normal distribution.
        /// </summary>
        /// <param name="x">The x value for the function.</param>
        /// <returns>The value f(x) at the specified position x.</returns>
        public static double NormalDistributionDensityFunction(double x)
        {
            return (0.5*Math.PI)*Math.Pow(Math.E, -0.5*(x*x));
        }

        /// <summary>
        /// The Zeta function.
        /// </summary>
        /// <param name="s">The exponent s for the zetafunction.</param>
        /// <returns>The value for Z(s).</returns>
        public static double ZetaFunction(double s)
        {
            return SpecialFunctions.ZetaFunction(s, SpecialFunctions.StandardIntegrationValue);
        }

        /// <summary>
        /// The Zeta function.
        /// </summary>
        /// <param name="s">The exponent s for the zetafunction.</param>
        /// <param name="limit">The limit of integration.</param>
        /// <returns>The value for Z(s).</returns>
        public static double ZetaFunction(double s, double limit)
        {
            double result = 1;
            double tempuri = 0;

            for (int i = 2; i < limit; i++)
            {
                tempuri = 1/(Math.Pow(i, s));
                result += tempuri;

                if (tempuri < limit)
                {
                    break;
                }
            }

            return result;
        }

        /// <summary>
        /// The Trigamma function.
        /// </summary>
        /// <param name="z">The z value for the function.</param>
        /// <returns>The value f(z) at the specified position z.</returns>
        public static double TrigammaFunction(double z)
        {
            return SpecialFunctions.TrigammaFunction(z, SpecialFunctions.StandardIntegrationValue);
        }

        /// <summary>
        /// The Trigamma function.
        /// </summary>
        /// <param name="z">The z value for the function.</param>
        /// <param name="limit">The limit of integration.</param>
        /// <returns>The value f(z) at the specified position z.</returns>
        public static double TrigammaFunction(double z, uint limit)
        {
            double result = 0;

            for (int n = 1; n <= limit; n++)
            {
                result += (1/Math.Pow(z + n, 2));
            }

            return result;
        }

        /// <summary>
        /// Computes a approximation result of the gamma function by the Stirling formula for 
        /// values x > 0.
        /// </summary>
        /// <param name="x">The specified number x, for which the gamma function should be calculated.</param>
        /// <returns>Returns the approximation value of the gamma function for the specified number. </returns>
        public static double GammaFunctionApproximation(int x)
        {
            return GammaFunctionApproximation((double) x);
        }

        /// <summary>
        /// Computes a approximation result of the gamma function by the Stirling formula for 
        /// values x > 0.
        /// </summary>
        /// <param name="x">The specified number x, for which the gamma function should be calculated.</param>
        /// <returns>Returns the approximation value of the gamma function for the specified number. </returns>
        public static double GammaFunctionApproximation(double x)
        {
            return Math.Sqrt(2*Math.PI)*Math.Pow(x, x - 0.5)*Math.Pow(Math.E, -x);
        }

        /// <summary>
        /// The Gamma function for real numbers.
        /// </summary>
        /// <param name="x">The value for evaluation.</param>
        /// <returns>The result of the real gamma function at position x.</returns>
        public static double GammaFunction(double x)
        {
            //Copyright (c) 2003, C. Bond. All rights reserved.
            int i, k, m;
            double ga = 0;
            double gr = 0;
            double r = 0;
            double z = 0;

            double[] g = {
                             1.0,
                             0.5772156649015329,
                             -0.6558780715202538,
                             -0.420026350340952e-1,
                             0.1665386113822915,
                             -0.421977345555443e-1,
                             -0.9621971527877e-2,
                             0.7218943246663e-2,
                             -0.11651675918591e-2,
                             -0.2152416741149e-3,
                             0.1280502823882e-3,
                             -0.201348547807e-4,
                             -0.12504934821e-5,
                             0.1133027232e-5,
                             -0.2056338417e-6,
                             0.6116095e-8,
                             0.50020075e-8,
                             -0.11812746e-8,
                             0.1043427e-9,
                             0.77823e-11,
                             -0.36968e-11,
                             0.51e-12,
                             -0.206e-13,
                             -0.54e-14,
                             0.14e-14
                         };

            if (x > 171.0)
            {
                return 1e308;
            }

            if (x == (int) x)
            {
                if (x > 0.0)
                {
                    ga = 1.0;
                    for (i = 2; i < x; i++)
                    {
                        ga *= i;
                    }
                }
                else
                {
                    ga = 1e308;
                }
            }
            else
            {
                if (ExtendedMath.Fabs(x) > 1.0)
                {
                    z = ExtendedMath.Fabs(x);
                    m = (int) z;
                    r = 1.0;

                    for (k = 1; k <= m; k++)
                    {
                        r *= (z - k);
                    }

                    z -= m;
                }
                else
                {
                    z = x;
                }

                gr = g[24];

                for (k = 23; k >= 0; k--)
                {
                    gr = gr*z + g[k];
                }

                ga = 1.0/(gr*z);

                if (ExtendedMath.Fabs(x) > 1.0)
                {
                    ga *= r;

                    if (x < 0.0)
                    {
                        ga = -Math.PI/(x*ga*Math.Sin(Math.PI*x));
                    }
                }
            }

            return ga;
        }

        /// <summary>
        /// The Gamma function for complex numbers.
        /// </summary>
        /// <param name="z">The value for evaluation.</param>
        /// <returns>The result of the complex gamma function at position z.</returns>
        public static ComplexNumber ComplexGammaFunction(ComplexNumber z)
        {
            //Based on the algorithm of C. Bond. Copyright (c) 2003, C. Bond. All rights reserved.
            double na = 0;
            double x1 = 0;
            double[] a =
                {
                    8.333333333333333e-02,
                    -2.777777777777778e-03,
                    7.936507936507937e-04,
                    -5.952380952380952e-04,
                    8.417508417508418e-04,
                    -1.917526917526918e-03,
                    6.410256410256410e-03,
                    -2.955065359477124e-02,
                    1.796443723688307e-01,
                    -1.39243221690590
                };

            double x = z.RealNumberPart;
            double y = z.ImaginaryNumberPart;

            if (x > 171)
            {
                return new ComplexNumber(1e308, 0);
            }
            if ((y == 0.0) && (x == (int) x) && (x <= 0.0))
            {
                return new ComplexNumber(1e308, 0);
            }
            if (x < 0.0)
            {
                x1 = x;
                x = -x;
                y = -y;
            }

            double x0 = x;

            if (x <= 7.0)
            {
                na = (int) (7.0 - x);
                x0 = x + na;
            }

            double q1 = Math.Sqrt(x0*x0 + y*y);
            double th = Math.Atan(y/x0);
            double gr = (x0 - 0.5)*Math.Log(q1) - th*y - x0 + 0.5*Math.Log(2.0*Math.PI);
            double gi = th*(x0 - 0.5) + y*Math.Log(q1) - y;

            for (int k = 0; k < 10; k++)
            {
                double t = Math.Pow(q1, -1.0 - 2.0*k);

                gr += (a[k]*t*Math.Cos((2.0*k + 1.0)*th));
                gi -= (a[k]*t*Math.Sign((2.0*k + 1.0)*th));
            }

            if (x <= 7.0)
            {
                double gr1 = 0.0;
                double gi1 = 0.0;

                for (int j = 0; j < na; j++)
                {
                    gr1 += (0.5*Math.Log((x + j)*(x + j) + y*y));
                    gi1 += Math.Atan(y/(x + j));
                }

                gr -= gr1;
                gi -= gi1;
            }

            if (x1 <= 0.0)
            {
                double th1 = Math.Atan(y/x);
                double sr = -Math.Sin(Math.PI*x)*Math.Cosh(Math.PI*y);
                double si = -Math.Cos(Math.PI*x)*Math.Sinh(Math.PI*y);
                double q2 = Math.Sqrt(sr*sr + si*si);
                double th2 = Math.Atan(si/sr);

                q1 = Math.Sqrt(x*x + y*y);

                if (sr < 0.0)
                {
                    th2 += Math.PI;
                }

                gr = Math.Log(Math.PI/(q1*q2)) - gr;
                gi = -th1 - th2 - gi;
            }

            double g0 = Math.Exp(gr);

            gr = g0*Math.Cos(gi);
            gi = g0*Math.Sin(gi);

            return new ComplexNumber(gr, gi);
        }

        /// <summary>
        /// The complex Loggamma function ln T(z).
        /// </summary>
        /// <param name="z">The value for evaluation.</param>
        /// <returns>The value ln T(z) at the specified position.</returns>
        public static ComplexNumber ComplexLogGammaFunction(ComplexNumber z)
        {
            //Based on the algorithm of C. Bond. Copyright (c) 2003, C. Bond. All rights reserved.
            double na = 0;
            double x1 = 0;
            double[] a =
                {
                    8.333333333333333e-02,
                    -2.777777777777778e-03,
                    7.936507936507937e-04,
                    -5.952380952380952e-04,
                    8.417508417508418e-04,
                    -1.917526917526918e-03,
                    6.410256410256410e-03,
                    -2.955065359477124e-02,
                    1.796443723688307e-01,
                    -1.39243221690590
                };

            double x = z.RealNumberPart;
            double y = z.ImaginaryNumberPart;

            if (x > 171)
            {
                return new ComplexNumber(1e308, 0);
            }
            if ((y == 0.0) && (x == (int) x) && (x <= 0.0))
            {
                return new ComplexNumber(1e308, 0);
            }
            if (x < 0.0)
            {
                x1 = x;
                x = -x;
                y = -y;
            }

            double x0 = x;

            if (x <= 7.0)
            {
                na = (int) (7.0 - x);
                x0 = x + na;
            }

            double q1 = Math.Sqrt(x0*x0 + y*y);
            double th = Math.Atan(y/x0);
            double gr = (x0 - 0.5)*Math.Log(q1) - th*y - x0 + 0.5*Math.Log(2.0*Math.PI);
            double gi = th*(x0 - 0.5) + y*Math.Log(q1) - y;

            for (int k = 0; k < 10; k++)
            {
                double t = Math.Pow(q1, -1.0 - 2.0*k);

                gr += (a[k]*t*Math.Cos((2.0*k + 1.0)*th));
                gi -= (a[k]*t*Math.Sign((2.0*k + 1.0)*th));
            }

            if (x <= 7.0)
            {
                double gr1 = 0.0;
                double gi1 = 0.0;

                for (int j = 0; j < na; j++)
                {
                    gr1 += (0.5*Math.Log((x + j)*(x + j) + y*y));
                    gi1 += Math.Atan(y/(x + j));
                }

                gr -= gr1;
                gi -= gi1;
            }

            if (x1 <= 0.0)
            {
                double th1 = Math.Atan(y/x);
                double sr = -Math.Sin(Math.PI*x)*Math.Cosh(Math.PI*y);
                double si = -Math.Cos(Math.PI*x)*Math.Sinh(Math.PI*y);
                double q2 = Math.Sqrt(sr*sr + si*si);
                double th2 = Math.Atan(si/sr);

                q1 = Math.Sqrt(x*x + y*y);

                if (sr < 0.0)
                {
                    th2 += Math.PI;
                }

                gr = Math.Log(Math.PI/(q1*q2)) - gr;
                gi = -th1 - th2 - gi;
            }

            return new ComplexNumber(gr, gi);
        }

        /// <summary>
        /// The lower incomplete Gammafunction T(a,x).
        /// </summary>
        /// <param name="a">The value for evaluation a.</param>
        /// <param name="x">The value for evaluation x.</param>
        /// <returns>
        /// The value ln T(a,x) at the specified position.
        /// </returns>
        public static double LowerIncompleteGammaFunction(double a, double x)
        {
            double result = 0;
            double tempuri1 = 0;
            double tempuri2 = 0;

            IncompleteGammaFunction(a, x, ref result, ref tempuri1, ref tempuri2);

            return result;
        }

        /// <summary>
        /// The upper incomplete Gammafunction T(a,x).
        /// </summary>
        /// <param name="a">The value for evaluation a.</param>
        /// <param name="x">The value for evaluation x.</param>
        /// <returns>
        /// The value ln T(a,x) at the specified position.
        /// </returns>
        public static double UpperIncompleteGammaFunction(double a, double x)
        {
            double result = 0;
            double tempuri1 = 0;
            double tempuri2 = 0;

            IncompleteGammaFunction(a, x, ref tempuri1, ref tempuri2, ref result);

            return result;
        }

        /// <summary>
        /// The incomplete Gammafunction T(a,x).
        /// </summary>
        /// <param name="a">The value for evaluation a.</param>
        /// <param name="x">The value for evaluation x.</param>
        /// <returns>
        /// The value ln T(a,x) at the specified position.
        /// </returns>
        public static double IncompleteGammaFunction(double a, double x)
        {
            double result = 0;

            if ((a < 0.0) || (x < 0))
            {
                return result;
            }

            double xam = -x + a*Math.Log(x);

            if ((xam > 700) || (a > 170.0))
            {
                return result;
            }

            if (x == 0.0)
            {
                return result;
            }

            if (x <= 1.0 + a)
            {
                double s = 1.0/a;
                double r = s;

                for (int k = 1; k <= 60; k++)
                {
                    r *= x/(a + k);
                    s += r;

                    if (Math.Abs(r/s) < 1e-15)
                    {
                        break;
                    }
                }

                result = (Math.Exp(xam)*s)/SpecialFunctions.GammaFunction(a);
            }
            else
            {
                double t0 = 0.0;

                for (int k = 60; k >= 1; k--)
                {
                    t0 = (k - a)/(1.0 + k/(x + t0));
                }

                result = 1.0 - (Math.Exp(xam)/(x + t0))/SpecialFunctions.GammaFunction(a);
            }

            return result;
        }

        /// <summary>
        /// The incomplete Gammafunction T(a,x).
        /// </summary>
        /// <param name="a">The value for evaluation a.</param>
        /// <param name="x">The value for evaluation x.</param>
        /// <param name="gin">Returns the value for the lower incomplete Gammafunction.</param>
        /// <param name="gim">Returns the value for the incomplete Gammafunction.</param>
        /// <param name="gip">Returns the value for the upper incomplete Gammafunction.</param>
        /// <returns>
        /// The value 0 if the computation was succeed otherwise the calculation failed.
        /// </returns>
        public static int IncompleteGammaFunction(double a, double x, ref double gin, ref double gim, ref double gip)
        {
            int k;
            double ga;

            if ((a < 0.0) || (x < 0))
            {
                return 1;
            }

            double xam = -x + a*Math.Log(x);

            if ((xam > 700) || (a > 170.0))
            {
                return 1;
            }

            if (x == 0.0)
            {
                gin = 0.0;
                gim = GammaFunction(a);
                gip = 0.0;
                return 0;
            }

            if (x <= 1.0 + a)
            {
                double s = 1.0/a;
                double r = s;

                for (k = 1; k <= 60; k++)
                {
                    r *= x/(a + k);
                    s += r;

                    if (Math.Abs(r/s) < 1e-15)
                    {
                        break;
                    }
                }

                gin = Math.Exp(xam)*s;
                ga = GammaFunction(a);
                gip = gin/ga;
                gim = ga - gin;
            }
            else
            {
                double t0 = 0.0;

                for (k = 60; k >= 1; k--)
                {
                    t0 = (k - a)/(1.0 + k/(x + t0));
                }

                gim = Math.Exp(xam)/(x + t0);
                ga = GammaFunction(a);
                gin = ga - gim;
                gip = 1.0 - gim/ga;
            }

            return 0;
        }

        /// <summary>
        /// The Loggamma function ln T(x).
        /// </summary>
        /// <param name="x">The value for evaluation.</param>
        /// <returns>The value ln T(x) at the specified position.</returns>
        public static double LogGammaFunction(double x)
        {
            //Based on the algorithm of C. Bond. Copyright (c) 2003, C. Bond. All rights reserved.
            int n = 0;
            double x0 = x;
            double[] a = new double[]
                             {
                                 8.333333333333333e-02,
                                 -2.777777777777778e-03,
                                 7.936507936507937e-04,
                                 -5.952380952380952e-04,
                                 8.417508417508418e-04,
                                 -1.917526917526918e-03,
                                 6.410256410256410e-03,
                                 -2.955065359477124e-02,
                                 1.796443723688307e-01,
                                 -1.39243221690590
                             };

            if (x <= 0.0)
            {
                return 1e308;
            }
            if ((x == 1.0) || (x == 2.0))
            {
                return 0.0;
            }
            if (x <= 7.0)
            {
                n = (int) (7 - x);
                x0 = x + n;
            }

            double x2 = 1.0/(x0*x0);
            double gl0 = a[9];
            const double xp = 2.0*Math.PI;

            for (int k = 8; k >= 0; k--)
            {
                gl0 = gl0*x2 + a[k];
            }

            double gl = gl0/x0 + 0.5*Math.Log(xp) + (x0 - 0.5)*Math.Log(x0) - x0;

            if (x <= 7.0)
            {
                for (int k = 1; k <= n; k++)
                {
                    gl -= Math.Log(x0 - 1.0);
                    x0 -= 1.0;
                }
            }

            return gl;
        }

        /// <summary>
        /// The Beta function.
        /// </summary>
        /// <param name="x">The first parameter x of the beta function.</param>
        /// <param name="y">The second parameter y of the beta function.</param>
        /// <returns>The result B(x,y) of the beta function.</returns>
        public static double BetaFunction(double x, double y)
        {
            return (SpecialFunctions.GammaFunction(x)*SpecialFunctions.GammaFunction(y))/
                   SpecialFunctions.GammaFunction(x + y);
        }

        /// <summary>
        /// The regularized incomplete beta function.
        /// </summary>
        /// <param name="x">The first paramter of the function.</param>
        /// <param name="a">The second paramter of the function.</param>
        /// <param name="b">The tird paramter of the function.</param>
        /// <returns>The value at the specified position x.</returns>
        public static double RegularizedIncompleteBetaFunction(double x, double a, double b)
        {
            if (x == 0.0 || x == 1.0)
            {
                return x;
            }
            else
            {
                bool recur;
                double f;

                if (x > 0.5)
                {
                    f = a;
                    a = b;
                    b = f;
                    x = 1.0 - x;
                    recur = true;
                }
                else
                {
                    recur = false;
                }

                int n = 1;
                int m = 0;
                double g = 0;
                double fn1 = 1;
                double fn2 = 0;
                double gn1 = 1;
                double gn2 = 1;
                double pq = a + b;
                bool neven = false;

                f = 1;

                do
                {
                    double dn;

                    if (neven)
                    {
                        m++;
                        dn = m*x*(b - m)/(a + n - 1.0)/(a + n);
                    }
                    else
                    {
                        dn = -x*(a + m)*(pq + m)/(a + n - 1.0)/(a + n);
                    }

                    g = f;

                    double fn = fn1 + dn*fn2;
                    double gn = gn1 + dn*gn2;

                    neven = (!neven);
                    f = fn/gn;
                    fn2 = fn1;
                    fn1 = fn;
                    gn2 = gn1;
                    gn1 = gn;
                    n++;
                } while (Math.Abs((f - g)/f) > 1e-10);

                f = f*Math.Pow(x, a)*Math.Pow(1.0 - x, b)*GammaFunction(a + b)/GammaFunction(a + 1.0)/GammaFunction(b);

                if (recur)
                {
                    f = 1.0 - f;
                }

                return f;
            }
        }

        /// <summary>
        /// The Central Beta function.
        /// </summary>
        /// <param name="x">The value for evaluation.</param>
        /// <returns>The result of Central Beta function at position x.</returns>
        public static double CentralBetaFunction(double x)
        {
            return SpecialFunctions.BetaFunction(x, x);
        }

        /// <summary>
        /// Computes a approximation result of the Beta function by the Stirling formula for 
        /// values x > 0.
        /// </summary>
        /// <param name="x">The first parameter x of the beta function.</param>
        /// <param name="y">The second parameter y of the beta function.</param>
        /// <returns>The approximated result B(x,y) of the beta function.</returns>
        public static double BetaFunctionApproximation(double x, double y)
        {
            return Math.Sqrt(2*Math.PI)*
                   ((Math.Pow(x, x - 0.5)*Math.Pow(y, y - 0.5))/(Math.Pow(x + y, x + y - 0.5)));
        }

        /// <summary>
        /// The Dirichlet beta function.
        /// </summary>
        /// <param name="s">The value for evaluation.</param>
        /// <param name="limit">The limit of integration.</param>
        /// <returns>The result of the Dirichlet beta function at position s.</returns>
        public static double DirichletBetaFunction(double s, int limit)
        {
            double result = 0;

            for (int n = 0; n < limit; n++)
            {
                result += (Math.Pow(-1, n))/(Math.Pow(2.0*n + 1.0, s));
            }

            return result;
        }

        /// <summary>
        /// The Hurwitz Zeta function, which is define for s > 1.
        /// </summary>
        /// <param name="s">The first parameter s of the function.</param>
        /// <param name="a">The second parameter a of the function.</param>
        /// <param name="limit">The limit of integration.</param>
        /// <returns>The result f(s,a) of the Hurwitz zeta function.</returns>
        public static double HurwitzZetaFunction(double s, double a, int limit)
        {
            double result = 0;

            for (int k = 0; k < limit; k++)
            {
                result += 1.0/Math.Pow(k + a, s);
            }

            return result;
        }

        /// <summary>
        /// The Lerch Transcendent.
        /// </summary>
        /// <param name="z">The first parameter z.</param>
        /// <param name="s">The second parameter s.</param>
        /// <param name="a">The third parameter a.</param>
        /// <param name="limit">The limit of integration.</param>
        /// <returns>The result f(z,s,a) of the Lerch Transcendent.</returns>
        public static double LerchTranscendent(double z, double s, double a, int limit)
        {
            double result = 0;

            for (int k = 0; k < limit; k++)
            {
                result += Math.Pow(z, k)/Math.Pow(k + a, s);
            }

            return result;
        }

        /// <summary>
        /// The Polygamma function.
        /// </summary>
        /// <param name="z">The first parameter z of the polygamma function.</param>
        /// <param name="n">The second parameter n of the polygamma function.</param>
        /// <param name="limit">The limit of integration.</param>
        /// <returns>The result fn(z) of the Polygamma function.</returns>
        public static double PolyGammaFunction(double z, uint n, int limit)
        {
            double result = 0;

            for (int k = 0; k < limit; k++)
            {
                result += 1.0/Math.Pow(z + k, n + 1);
            }

            return Math.Pow(-1, n + 1)*Factorials.Factorial(n)*result;
        }

        /// <summary>
        /// The Totient function
        /// </summary>
        /// <param name="x">The value for evaluation.</param>
        /// <returns>The result of the function for the specified position x.</returns>
        public static uint TotientFunction(int x)
        {
            uint result = 0;

            for (int i = 0; i < x; i++)
            {
                result += (uint) Math.Floor(1.0/ExtendedMath.GCD(i, x));
            }

            return result;
        }

        /// <summary>
        /// The Totient summatory function.
        /// </summary>
        /// <param name="x">The specified value x.</param>
        /// <returns>The sum of all Totient function up to the specified value x.</returns>
        public static UInt64 TotientSummatoryFunction(int x)
        {
            UInt64 result = 0;

            for (int i = 0; i <= x; i++)
            {
                result += SpecialFunctions.TotientFunction(i);
            }

            return result;
        }

        /// <summary>
        /// The divisor function.
        /// </summary>
        /// <param name="n">The specified value n for evaluation.</param>
        /// <param name="x">The xth powers of the divisors of n.</param>
        /// <returns>The result of the divisor function ox(n).</returns>
        public static uint DivisorFunction(int n, int x)
        {
            double result = 0;

            for (uint u = 1; u <= n; u++)
            {
                double tempuri = 0;

                for (int v = 1; v <= u; v++)
                {
                    tempuri += Math.Cos((2*Math.PI*v*n)/u);
                }

                result += Math.Pow(u, x - 1)*tempuri;
            }

            return (uint) Math.Ceiling(result);
        }

        /// <summary>
        /// The divisor function.
        /// </summary>
        /// <param name="n">The specified value n for evaluation.</param>
        /// <param name="x">The xth powers of the divisors of n.</param>
        /// <returns>The result of the divisor function ox(n).</returns>
        public static Int64 DivisorFunction(Int64 n, Int64 x)
        {
            Int64 result = 0;

            for (uint u = 1; u <= n; u++)
            {
                Int64 tempuri = 0;

                for (int v = 1; v <= u; v++)
                {
                    tempuri += (Int64) Math.Cos((2*Math.PI*v*n)/u);
                }

                result += (Int64) Math.Pow(u, x - 1)*tempuri;
            }

            return result;
        }

        /// <summary>
        /// The Sigma function.
        /// </summary>
        /// <param name="n">The specified value n for evaluation.</param>
        /// <returns>The result of the Sigma function o1(n).</returns>
        public static uint SigmaFunction(int n)
        {
            return SpecialFunctions.DivisorFunction(n, 1);
        }

        /// <summary>
        ///Approximation of the sine function.
        /// </summary>
        /// <param name="x">The value for evaluation.</param>
        /// <returns>The result of the approximation.</returns>
        public static double SineApproximation(double x)
        {
            return x - (Math.Pow(x, 3)/6) + (Math.Pow(x, 5)/120) - (Math.Pow(x, 7)/5040);
        }

        /// <summary>
        /// Approximation of the sine function.
        /// </summary>
        /// <param name="x">The value for evaluation.</param>
        /// <param name="iterations">The number of iterations for the approximation.</param>
        /// <returns>The result of the approximation.</returns>
        public static double SineApproximation(double x, int iterations)
        {
            double result = 0;
            UInt64 end = (UInt64) iterations;

            for (UInt64 i = 0; i <= end; i++)
            {
                result += (Math.Pow(-1, i)/Factorials.Factorial(2*i + 1))*Math.Pow(x, 2*i + 1);
            }

            return result;
        }

        /// <summary>
        /// Cosines the approximation.
        /// </summary>
        /// <param name="x">The value for evaluation.</param>
        /// <param name="iterations">The number of iterations for the approximation.</param>
        /// <returns>The result of the approximation.</returns>
        public static double CosineApproximation(double x, int iterations)
        {
            UInt64 end = (UInt64) iterations;
            double result = 0;

            for (UInt64 i = 0; i <= end; i++)
            {
                result += (Math.Pow(-1, i)/Factorials.Factorial(2*i))*Math.Pow(x, 2*i);
            }

            return result;
        }

        /// <summary>
        /// The Errorfunction. 
        /// </summary>
        /// <param name="z">The value for evaluation.</param>
        /// <returns>The result of the function at the specified position z.</returns>
        public static double ErrorFunction(double z)
        {
            double t = 1.0/(1.0 + 0.47047*Math.Abs(z));
            double poly = t*(0.3480242 + t*(-0.0958798 + t*(0.7478556)));
            double ans = 1.0 - poly*Math.Exp(-z*z);

            if (z >= 0)
            {
                return ans;
            }

            return -ans;
        }

        /// <summary>
        /// The complementary Errorfunction.
        /// </summary>
        /// <param name="z">The value for evaluation.</param>
        /// <returns>The result of the function at the specified position z.</returns>
        public static double ComplementaryErrorFunction(double z)
        {
            return 1 - SpecialFunctions.ErrorFunction(z);
        }

        /// <summary>
        /// The WhittakerM function Mk,m(z). Its the Whittaker function of the first kind.
        /// </summary>
        /// <param name="k">The value of k.</param>
        /// <param name="m">The value of m.</param>
        /// <param name="z">The value of z.</param>
        /// <returns>The result of the function at the specified position.</returns>
        public static double WhittakerMFunction(int k, int m, double z)
        {
            double result = 0;

            for (uint n = 0; n <= 10; n++)
            {
                result += (SpecialFunctions.PochhammerSymbol(m - k + 0.5, n)
                           /(Factorials.Factorial(n)*SpecialFunctions.PochhammerSymbol(2*m + 1, n)))
                          *Math.Pow(z, n);
            }

            return Math.Pow(z, 0.5 + m)*Math.Pow(Math.E, -z/2)*result;
        }

        /// <summary>
        /// The confluent hypergeometric function F1(a; b; z) of the first kind.
        /// </summary>
        /// <param name="a">The value of a.</param>
        /// <param name="b">The value of b.</param>
        /// <param name="z">The value of z.</param>
        /// <returns>The result of the function at the specified position.</returns>
        public static double ConfluentHypergeometricFunctionFirstKind(double a, double b, double z)
        {
            double result = 0;
            UInt64[] factorials = Factorials.FactorialList20();

            for (UInt64 k = 0; k < 20; k++)
            {
                result += (SpecialFunctions.PochhammerSymbol(a, k)/SpecialFunctions.PochhammerSymbol(b, k))
                          *(Math.Pow(z, k)/factorials[k]);
            }

            return result;
        }

        /// <summary>
        /// The generalized hypergeometric function.
        /// </summary>
        /// <param name="a">The first series a.</param>
        /// <param name="b">The second series b.</param>
        /// <param name="z">The value for evaluation.</param>
        /// <param name="limit">The limit of integration, but maximum 65 to prevent performance leaks.</param>
        /// <returns>The result of the function at the specified position z.</returns>
        public static double GeneralizedHypergeometricFunction(double[] a, double[] b, double z, int limit)
        {
            if (limit > 65)
            {
                limit = 65;
            }

            double result = 0;
            double resultA = 1;
            double resultB = 1;
            UInt64[] factorials = Factorials.FactorialList20();

            if (limit > 20)
            {
                factorials = Factorials.FactorialList65();
            }

            for (int k = 0; k < limit; k++)
            {
                resultA = resultB = 1;

                for (int i = 0; i < a.Length; i++)
                {
                    resultA *= SpecialFunctions.PochhammerSymbol(a[i], k);
                }

                for (int i = 0; i < b.Length; i++)
                {
                    resultB *= SpecialFunctions.PochhammerSymbol(b[i], k);
                }

                result += (resultA/resultB)*(Math.Pow(z, k)/factorials[k]);
            }

            return result;
        }

        /// <summary>
        /// The Pochhammer symbol (x)n.
        /// </summary>
        /// <param name="x">The value for evaluation.</param>
        /// <param name="n">The order n for the Pochhammer symbol.</param>
        /// <returns>The result of the function at the specified position x.</returns>
        public static double PochhammerSymbol(double x, double n)
        {
            return SpecialFunctions.GammaFunction(x + n)/SpecialFunctions.GammaFunction(x);
        }

        /// <summary>
        /// The Kummerfunction 2F1(1/2 + m - k, -n; 2m + 1; 1) of the first kind.
        /// </summary>
        /// <param name="k">The value of k.</param>
        /// <param name="m">The value of m.</param>
        /// <param name="n">The value of n.</param>
        /// <returns>The result of the function at the specified position.</returns>
        public static double KummerFunctionFirstKind(double k, double m, double n)
        {
            double a = SpecialFunctions.GammaFunction(2*m + 1);
            double b = SpecialFunctions.GammaFunction(m + 0.5 + k + n);
            double c = SpecialFunctions.GammaFunction(m + 0.5 + k);
            double d = SpecialFunctions.GammaFunction(2*m + 1 + n);

            return (a*b)/(c*d);
        }

        /// <summary>
        /// Bessel function Jn(x) of the first kind.
        /// </summary>
        /// <param name="x">The value for evaluation.</param>
        /// <param name="n">The order n of the function.</param>
        /// <returns>The result of the function at the specified position x.</returns>
        public static double BesselFunctionFirstKind(double x, int n)
        {
            bool isNull = false;

            if (n == 0)
            {
                n++;
                isNull = true;
            }

            double[] jn = new double[n + 1];
            double[] yn = new double[n + 1];
            double[] jnp = new double[n + 1];
            double[] ynp = new double[n + 1];

            SpecialFunctions.BesselFunction(n, x, jn, yn, jnp, ynp);

            if (isNull)
            {
                return jn[0];
            }

            return jn[n];
        }

        /// <summary>
        /// Bessel function Yn(x) of the second kind.
        /// </summary>
        /// <param name="x">The value for evaluation.</param>
        /// <param name="n">The order n of the function.</param>
        /// <returns>The result of the function at the specified position x.</returns>
        public static double BesselFunctionSecondKind(double x, int n)
        {
            bool isNull = false;

            if (n==0)
            {
                n++;
                isNull = true;
            }

            double[] jn = new double[n + 1];
            double[] yn = new double[n + 1];
            double[] jnp = new double[n + 1];
            double[] ynp = new double[n + 1];

            SpecialFunctions.BesselFunction(n, x, jn, yn, jnp, ynp);

            if (isNull)
            {
                return yn[0];
            }

            return yn[n];
        }

        /// <summary>
        /// Computes the Bessel function of the first and second kind.
        /// </summary>
        /// <param name="x">The parameter of the Bessel function.</param>
        /// <param name="j0">Bessel function of 1st kind, 0th order.</param>
        /// <param name="j1">Bessel function of 1st kind, 1st order.</param>
        /// <param name="y0">Bessel function of 2nd kind, 0th order.</param>
        /// <param name="y1">Bessel function of 2nd kind, 1st order.</param>
        /// <param name="j0p">Derivative of Bessel function of 1st kind, 0th order.</param>
        /// <param name="j1p">Derivative of Bessel function of 1st kind, 1st order.</param>
        /// <param name="y0p">Derivative of Bessel function of 2nd kind, 0th order.</param>
        /// <param name="y1p">Derivative of Bessel function of 2nd kind, 1st order.</param>
        /// <returns>Returns 0, if the computation of the Bessel function was correct otherwise, false.</returns>
        public static int BesselFunction(double x, ref double j0, ref double j1, ref double y0, ref double y1, ref double j0p,
                                    ref double j1p, ref double y0p, ref double y1p)
        {
            double M_2_PI = 2 / Math.PI;
            double M_PI_2 = Math.PI / 2;
            double M_PI_4 = Math.PI / 4;
            double eps = 1e-15;
            double el = 0.5772156649015329;
            int k;
            double[] a = {
                             -7.03125e-2,
                             0.112152099609375,
                             -0.5725014209747314,
                             6.074042001273483,
                             -1.100171402692467e2,
                             3.038090510922384e3,
                             -1.188384262567832e5,
                             6.252951493434797e6,
                             -4.259392165047669e8,
                             3.646840080706556e10,
                             -3.833534661393944e12,
                             4.854014686852901e14,
                             -7.286857349377656e16,
                             1.279721941975975e19
                         };

            double[] b = {
                             7.32421875e-2,
                             -0.2271080017089844,
                             1.727727502584457,
                             -2.438052969955606e1,
                             5.513358961220206e2,
                             -1.825775547429318e4,
                             8.328593040162893e5,
                             -5.006958953198893e7,
                             3.836255180230433e9,
                             -3.649010818849833e11,
                             4.218971570284096e13,
                             -5.827244631566907e15,
                             9.476288099260110e17,
                             -1.792162323051699e20
                         };

            double[] a1 = {
                              0.1171875,
                              -0.1441955566406250,
                              0.6765925884246826,
                              -6.883914268109947,
                              1.215978918765359e2,
                              -3.302272294480852e3,
                              1.276412726461746e5,
                              -6.656367718817688e6,
                              4.502786003050393e8,
                              -3.833857520742790e10,
                              4.011838599133198e12,
                              -5.060568503314727e14,
                              7.572616461117958e16,
                              -1.326257285320556e19
                          };

            double[] b1 = {
                              -0.1025390625,
                              0.2775764465332031,
                              -1.993531733751297,
                              2.724882731126854e1,
                              -6.038440767050702e2,
                              1.971837591223663e4,
                              -8.902978767070678e5,
                              5.310411010968522e7,
                              -4.043620325107754e9,
                              3.827011346598605e11,
                              -4.406481417852278e13,
                              6.065091351222699e15,
                              -9.833883876590679e17,
                              1.855045211579828e20
                          };

            if (x < 0.0)
            {
                return 1;
            }

            if (x == 0.0)
            {
                j0 = 1.0;
                j1 = 0.0;
                y0 = -1e308;
                y1 = -1e308;
                j0p = 0.0;
                j1p = 0.5;
                y0p = 1e308;
                y1p = 1e308;

                return 0;
            }

            double x2 = x * x;

            if (x <= 12.0)
            {
                j0 = 1.0;

                double r = 1.0;

                for (k = 1; k <= 30; k++)
                {
                    r *= -0.25 * x2 / (k * k);
                    j0 += r;

                    if (Math.Abs(r) < Math.Abs(j0) * 1e-15)
                    {
                        break;
                    }
                }

                j1 = 1.0;
                r = 1.0;

                for (k = 1; k <= 30; k++)
                {
                    r *= -0.25 * x2 / (k * (k + 1));
                    j1 += r;

                    if (Math.Abs(r) < Math.Abs(j1) * 1e-15)
                    {
                        break;
                    }
                }

                j1 *= 0.5 * x;

                double ec = Math.Log(0.5 * x) + el;
                double cs0 = 0.0;
                double w0 = 0.0;
                double r0 = 1.0;

                for (k = 1; k <= 30; k++)
                {
                    w0 += 1.0 / k;
                    r0 *= -0.25 * x2 / (k * k);
                    r = r0 * w0;
                    cs0 += r;

                    if (Math.Abs(r) < Math.Abs(cs0) * 1e-15)
                    {
                        break;
                    }
                }

                y0 = M_2_PI * (ec * j0 - cs0);

                double cs1 = 1.0;
                double w1 = 0.0;
                double r1 = 1.0;

                for (k = 1; k <= 30; k++)
                {
                    w1 += 1.0 / k;
                    r1 *= -0.25 * x2 / (k * (k + 1));
                    r = r1 * (2.0 * w1 + 1.0 / (k + 1));
                    cs1 += r;

                    if (Math.Abs(r) < Math.Abs(cs1) * 1e-15)
                    {
                        break;
                    }
                }

                y1 = M_2_PI * (ec * j1 - 1.0 / x - 0.25 * x * cs1);
            }
            else
            {
                int kz;

                if (x >= 50.0)
                {
                    kz = 8;
                }
                else if (x >= 35.0)
                {
                    kz = 10;
                }
                else
                {
                    kz = 12;
                }

                double t1 = x - M_PI_4;
                double p0 = 1.0;
                double q0 = -0.125 / x;

                for (k = 0; k < kz; k++)
                {
                    p0 += a[k] * Math.Pow(x, -2 * k - 2);
                    q0 += b[k] * Math.Pow(x, -2 * k - 3);
                }

                double cu = Math.Sqrt(M_2_PI / x);

                j0 = cu * (p0 * Math.Cos(t1) - q0 * Math.Sin(t1));
                y0 = cu * (p0 * Math.Sin(t1) + q0 * Math.Cos(t1));

                double t2 = x - 0.75 * Math.PI;
                double p1 = 1.0;
                double q1 = 0.375 / x;

                for (k = 0; k < kz; k++)
                {
                    p1 += a1[k] * Math.Pow(x, -2 * k - 2);
                    q1 += b1[k] * Math.Pow(x, -2 * k - 3);
                }

                j1 = cu * (p1 * Math.Cos(t2) - q1 * Math.Sin(t2));
                y1 = cu * (p1 * Math.Sin(t2) + q1 * Math.Cos(t2));
            }

            j0p = -j1;
            j1p = j0 - j1 / x;
            y0p = -y1;
            y1p = y0 - y1 / x;

            return 0;
        }

        /// <summary>
        /// Computes the Bessel function of the first and second kind.
        /// </summary>
        /// <param name="n">The order of the Bessel function.</param>
        /// <param name="x">Argument of Bessel function of 1st and 2nd kind.</param>
        /// <param name="jn">Bessel function of 1st kind, orders from 0 to nm.</param>
        /// <param name="yn">Bessel function of 2nd kind, orders from 0 to nm.</param>
        /// <param name="jnp">Derivative of Bessel function of 1st kind, orders from 0 to nm.</param>
        /// <param name="ynp">Derivative of Bessel function of 2nd kind, orders from 0 to nm.</param>
        /// <returns>Returns 0, if the computation of the Bessel function was correct otherwise, false.</returns>
        public static int BesselFunction(int n, double x, double[] jn, double[] yn, double[] jnp, double[] ynp)
        {
            double f = 1;
            double f1;
            int k;
            int nm;
            nm = n;

            if ((x < 0.0) || (n < 0))
            {
                return 1;
            }

            if (x < 1e-15)
            {
                int i;

                for (i = 0; i <= n; i++)
                {
                    jn[i] = 0.0;
                    yn[i] = -1e308;
                    jnp[i] = 0.0;
                    ynp[i] = 1e308;
                }

                jn[0] = 1.0;
                jnp[1] = 0.5;

                return 0;
            }

            SpecialFunctions.BesselFunction(x, ref jn[0], ref jn[1], ref yn[0], ref yn[1], ref jnp[0], ref jnp[1], ref ynp[0], ref ynp[1]);

            if (n < 2)
            {
                return 0;
            }

            double bj0 = jn[0];
            double bj1 = jn[1];

            if (n < (int)0.9 * x)
            {
                for (k = 2; k <= n; k++)
                {
                    jn[k] = 2.0 * (k - 1.0) * bj1 / x - bj0;
                    bj0 = bj1;
                    bj1 = jn[k];
                }
            }
            else
            {
                int m = BesselFunctionHelper1(x, 200);

                if (m < n)
                {
                    nm = m;
                }
                else
                {
                    m = BesselFunctionHelper2(x, n, 15);
                }

                double f2 = 0.0;

                f1 = 1.0e-100;

                for (k = m; k >= 0; k--)
                {
                    f = 2.0 * (k + 1.0) / x * f1 - f2;

                    if (k <= nm)
                    {
                        jn[k] = f;
                    }

                    f2 = f1;
                    f1 = f;
                }

                double cs;

                if (Math.Abs(bj0) > Math.Abs(bj1))
                {
                    cs = bj0 / f;
                }
                else
                {
                    cs = bj1 / f2;
                }
                for (k = 0; k <= nm; k++)
                {
                    jn[k] *= cs;
                }
            }

            for (k = 2; k <= nm; k++)
            {
                jnp[k] = jn[k - 1] - k * jn[k] / x;
            }

            double f0 = yn[0];

            f1 = yn[1];

            for (k = 2; k <= nm; k++)
            {
                f = 2.0 * (k - 1.0) * f1 / x - f0;
                yn[k] = f;
                f0 = f1;
                f1 = f;
            }

            for (k = 2; k <= nm; k++)
            {
                ynp[k] = yn[k - 1] - k * yn[k] / x;
            }

            return 0;
        }

        internal static int BesselFunctionHelper1(double x, int mp)
        {
            int i;
            int nn = 0;
            double a0 = Math.Abs(x);
            int n0 = (int)(1.1 * a0) + 1;
            double f0 = 0.5 * Math.Log10(6.28 * n0) - n0 * Math.Log10(1.36 * a0 / n0) - mp;
            int n1 = n0 + 5;
            double f1 = 0.5 * Math.Log10(6.28 * n1) - n1 * Math.Log10(1.36 * a0 / n1) - mp;

            for (i = 0; i < 20; i++)
            {
                nn = (int)(n1 - (n1 - n0) / (1.0 - f0 / f1));

                double f = 0.5 * Math.Log10(6.28 * nn) - nn * Math.Log10(1.36 * a0 / nn) - mp;

                if (Math.Abs(nn - n1) < 1)
                {
                    break;
                }
                n0 = n1;
                f0 = f1;
                n1 = nn;
                f1 = f;
            }

            return nn;
        }

        internal static int BesselFunctionHelper2(double x, int n, int mp)
        {
            double obj;
            int i, n0, nn = 0;
            double a0 = Math.Abs(x);
            double hmp = 0.5 * mp;
            double ejn = 0.5 * Math.Log10(6.28 * n) - n * Math.Log10(1.36 * a0 / n);

            if (ejn <= hmp)
            {
                obj = mp;
                n0 = (int)(1.1 * a0);

                if (n0 < 1)
                {
                    n0 = 1;
                }
            }
            else
            {
                obj = hmp + ejn;
                n0 = n;
            }

            double f0 = 0.5 * Math.Log10(6.28 * n0) - n0 * Math.Log10(1.36 * a0 / n0) - obj;
            int n1 = n0 + 5;
            double f1 = 0.5 * Math.Log10(6.28 * n1) - n1 * Math.Log10(1.36 * a0 / n1) - obj;

            for (i = 0; i < 20; i++)
            {
                nn = (int)(n1 - (n1 - n0) / (1.0 - f0 / f1));

                double f = 0.5 * Math.Log10(6.28 * nn) - nn * Math.Log10(1.36 * a0 / nn) - obj;

                if (Math.Abs(nn - n1) < 1)
                {
                    break;
                }

                n0 = n1;
                f0 = f1;
                n1 = nn;
                f1 = f;
            }

            return nn + 10;
        }

        /// <summary>
        /// The PSI function.
        /// </summary>
        /// <param name="x">The value for evaluation.</param>
        /// <returns>The result of the function at the specified position x.</returns>
        public static double PSIFunction(double x)
        {
            int n;
            double ps, x2;

            double[] a = {
                             -0.8333333333333e-01,
                             0.83333333333333333e-02,
                             -0.39682539682539683e-02,
                             0.41666666666666667e-02,
                             -0.75757575757575758e-02,
                             0.21092796092796093e-01,
                             -0.83333333333333333e-01,
                             0.4432598039215686
                         };

            double xa = Math.Abs(x);
            double s = 0.0;
            const double el = 0.5772156649015329;

            if ((x == (int) x) && (x <= 0.0))
            {
                ps = 1e308;

                return ps;
            }

            if (xa == (int) xa)
            {
                n = (int) xa;

                for (int k = 1; k < n; k++)
                {
                    s += 1.0/k;
                }

                ps = s - el;
            }
            else if ((xa + 0.5) == ((int) (xa + 0.5)))
            {
                n = (int) (xa - 0.5);

                for (int k = 1; k <= n; k++)
                {
                    s += 1.0/(2.0*k - 1.0);
                }

                ps = 2.0*s - el - 1.386294361119891;
            }
            else
            {
                if (xa < 10.0)
                {
                    n = 10 - (int) xa;

                    for (int k = 0; k < n; k++)
                    {
                        s += 1.0/(xa + k);
                    }

                    xa += n;
                }

                x2 = 1.0/(xa*xa);
                ps = Math.Log(xa) - 0.5/xa +
                     x2*
                     (((((((a[7]*x2 + a[6])*x2 + a[5])*x2 + a[4])*x2 + a[3])*x2 + a[2])*x2 + a[1])*x2 +
                      a[0]);
                ps -= s;
            }

            if (x < 0.0)
            {
                ps -= Math.PI*Math.Cos(Math.PI*x)/Math.Sin(Math.PI*x) - 1.0/x;
            }

            return ps;
        }

        /// <summary>
        /// The complex PSI function.
        /// </summary>
        /// <param name="z">The value for evaluation.</param>
        /// <returns>The result of the function at the specified position z.</returns>
        public static ComplexNumber ComplexPSIFunction(ComplexNumber z)
        {
            double y1 = 0;
            double th = 0;
            ComplexNumber ps;

            double[] a = {
                             -0.8333333333333e-01,
                             0.83333333333333333e-02,
                             -0.39682539682539683e-02,
                             0.41666666666666667e-02,
                             -0.75757575757575758e-02,
                             0.21092796092796093e-01,
                             -0.83333333333333333e-01,
                             0.4432598039215686
                         };

            double x = z.RealNumberPart;
            double y = z.ImaginaryNumberPart;
            double x1 = 0.0;
            int n = 0;

            if ((y == 0.0) && (x == (int) x) && (x <= 0.0))
            {
                ps = new ComplexNumber(1e308, 0);
            }
            else
            {
                if (x < 0.0)
                {
                    x1 = x;
                    y1 = y;
                    x = -x;
                    y = -y;
                }

                double x0 = x;

                if (x < 8.0)
                {
                    n = 8 - (int) x;
                    x0 = x + n;
                }

                if ((x0 == 0.0) && (y != 0.0))
                {
                    th = 0.5*Math.PI;
                }

                if (x0 != 0.0)
                {
                    th = Math.Atan(y/x0);
                }

                double q2 = x0*x0 + y*y;
                double q0 = Math.Sqrt(q2);
                double psr = Math.Log(q0) - 0.5*x0/q2;
                double psi = th + 0.5*y/q2;

                for (int k = 1; k <= 8; k++)
                {
                    psr += (a[k - 1]*Math.Pow(q2, -k)*Math.Cos(2.0*k*th));
                    psi -= (a[k - 1]*Math.Pow(q2, -k)*Math.Sin(2.0*k*th));
                }

                if (x < 8.0)
                {
                    double rr = 0.0;
                    double ri = 0.0;

                    for (int k = 1; k <= n; k++)
                    {
                        rr += ((x0 - k)/(Math.Pow(x0 - k, 2.0) + y*y));
                        ri += (y/(Math.Pow(x0 - k, 2.0) + y*y));
                    }

                    psr -= rr;
                    psi += ri;
                }

                if (x1 < 0.0)
                {
                    double tn = Math.Tan(Math.PI*x);
                    double tm = Math.Tan(Math.PI*y);
                    double ct2 = tn*tn + tm*tm;

                    psr = psr + x/(x*x + y*y) + Math.PI*(tn - tn*tm*tm)/ct2;
                    psi = psi - y/(x*x + y*y) - Math.PI*tm*(1.0 + tn*tn)/ct2;
                    x = x1;
                    y = y1;
                }

                ps = new ComplexNumber(psr, psi);
            }

            return ps;
        }
    }
}